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ABSTRACT 

Integrals  o£  a  function  of  a  single  variable  can  be  eiqnressed  as  the  sum  of 
a  numerical  quadrature  rule  and  a  remainder  term.  The  quadrature  rule  is  a 
linear  combination  of  function  values  and  weights,  or  the  integral  of  a  Taylor 
polynomial,  while  the  remainder  term  depends  on  some  derivative  of  the  integrand 
evaluated  at  an  unknown  point  in  the  interval  of  integration.  Numerical 
quadrature  is  made  self -validating  by  using  interval  computation  to  capture  both 
the  roundoff  and  truncation  errors  made  when  using  a  given  rule.  Necessary 
derivatives  can  be  generated  automatically  by  using  well-known  recurrence 
relations  for  Taylor  coefficients.  In  order  for  quadrature  methods  of  this  type 
to  be  accurate^ (in  the  sense  that  small  intervals  containing  the  exact  result 
are  produced)^  and  efficient  (to  obtain  results  of  given  accuracy  in  a  reasonably 
short  time  f',  an  accurate  scalar  product  and  an  adaptive  strategy  are  required. 
The  necessary  scalar  product  and  support  for  interval  arithmetic  are  provided  in 
Pascal-SC  (for  microcomputers )  and  ACRITH  (for  IBM  370  computers).  The  adaptive 
strategy  chooses  the  subintervals  of  integration  and  the  order  of  the  quadrature 
formula  used  in  each  subinterval  on  the  basis  of  guaranteed,  rather  than 
estimated,  information  about  the  error  of  the  numerical  integration  in  each 
subinterval.  The  program  described  'in  this  report-' implements  standard  Newton- 
Cotes,  Gaussian,  and  Taylor  series  methods  for  numerical  integration.  Ways  to 
handle  singularities  are  discussed,  and  comparisons  sure  given  with  a  standard 
numerical  integration  method. 
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SIGNIFICANCE  AND  EXPLANATION 


Routines  for  numerical  integration  are  among  the  most  heavily  used  programs 
in  computer  libraries  supporting  scientific,  engineering,  and  statistical 
computation.  The  ones  in  common  use  at  the  present  time,  such  as  CADRE  and 
QUADPACK,  produce  approximations  to  integrals  with  only  estimates  of  the  error 
in  the  value  returned.  The  program  described  in  this  report  performs  self¬ 
validating  numerical  quadrature,  returning  an  interval  in  which  the  desired 
integral  is  guaranteed  to  lie.  To  do  this,  automatic  differentiation  and 
interval  computation  are  used  to  capture  both  the  roundoff  and  truncation  error 
inherent  in  standard  Newton-Cotes,  Gaussian,  and  Taylor  polynomial  methods  for 
numerical  integration.  The  actual  integration  can  be  expressed  as  the  interval 
scalar  product  of  a  vector  of  function  values  with  a  vector  of  weights,  so  in 
order  to  obtain  accuracy  (small  intervals)  the  confutation  has  to  be  supported 
by  an  accurate  scalar  product  as  well  as  interval  arithmetic.  Fortunately, 
these  capabilities  are  provided  in  Pascal-SC  for  microcomputers,  and  ACRITH  for 
IBM  370  mainframe  coufuters.  The  program  described  in  this  report  is  mitten  in 
FORTRAN,  and  uses  ACRITH.  This  avoids  the  kind  of  special  inf  lementation  of 
interval  confutation  which  was  necessary  on  outdated  systems.  Since  the  IEEE 
standard  for  floating-point  arithmetic  requires  the  support  of  interval 
arithmetic,  it  can  be  assumed  that  self-validating  confutation  in  general  will 
be  much  easier  to  perform  as  conforming  systems  become  available.  Confarisons 
of  the  program  described  in  this  paper  with  the  standard  integration  routine 
QUADPACK  indicate  comparable  execution  times;  however,  the  results  of  QUADPACK 
lack  any  guarantee  of  accuracy. 
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ADAPTIVE.  SELF-VALIDATING  NUMERICAL  QUADRATURE 


George  F.  Corliss  and  L.  B.  Rail 


1.  Requirements  for  Automatic  Integration  Algorithms 

In  [2],  de  Boor  formulates  fundamental  requirements  for  an  automatic  algorithm  for  nu¬ 
merical  approximation  of  the  integral 


(1.1) 


of  a  function  of  a  single  real  variable.  Such  an  algorithm  requires  (i)  the  limits  of  integration 
a,  6,  (ii)  access  to  a  procedure  for  the  evaluation  of  /( x)  for  x  in  the  interval  of  integration, 
(iii)  tolerances  a.,p  on  the  desired  absolute  and  relative  error,  respectively,  and  (iv)  a  limit 
M  on  the  number  of  function  evaluations  allowed. 

As  output,  the  algorithm  should  produce  an  estimate  I*  for  the  value  of  If  which 
satisfies 


(1.2) 


I //  -  /*!  <  max{a,/>|//|}. 


Furthermore,  the  algorithm  should  be  efficient,  computing  as  few  function  values  as  pos¬ 
sible.  It  should  also  be  reliable,  which  will  be  taken  here  to  mean  that  either  the  desired 
accuracy  (1.2)  is  guaranteed,  or  a  message  to  the  contrary  is  returned  to  the  user,  possibly 
with  additional  information  about  the  cause  of  failure.  As  pointed  out  by  de  Boor  [2], 
algorithms  which  use  only  values  of  f(x)  at  a  finite  number  of  points  cannot  meet  the 
above  requirements  in  general;  nevertheless,  accurate  and  efficient  automatic  integration 
algorithms  can  be  formulated  for  wide  classes  of  integrands  |3j,  [231. 

This  paper  presents  automatic  quadrature  algorithms  which  attain  the  goals  of  reli¬ 
ability  and  efficiency  by  use  of  automatic  differentiation  and  interval  computation.  They 

Sponsored  in  part  by  the  U.  S.  Army  under  Contract  No.  DAAG29-80-C-0041. 


V  /  «*. 


make  use  of  information  about  the  integrand  on  entire  subintervals  of  integration,  rather 
than  at  a  discrete  set  of  points.  The  results  combine  the  self-validating  algorithms  of  Gray 
and  Rail  [8],  (9j,  [10],  and  the  notion  of  adaptive  quadrature  [2j,  [3],  [23].  Adaptation 
is  carried  out  on  the  basis  of  guaranteed,  rather  than  estimated,  bounds  for  the  error  of 
the  approximate  integration  over  each  subintervai.  Furthermore,  the  given  algorithm  has 
the  ability  to  detect  and  handle  certain  types  of  singularities  in  the  integrand,  and  even 
to  verify  nonexistence  of  the  integral  in  some  cases.  In  the  terminology  of  Rice  [23],  the 
method  described  here  has  the  following  features: 

Interval  Processor  Component: 

Variable  order  rule  with  remainder  using  interval 

arithmetic  to  give  guaranteed  bounds. 

Bound  Estimator  Component: 

Direct  analysis. 

• 

Special  Behavior  Components: 

Polynomials. 

Roundoff  level. 

Singularities  in  derivatives. 

Jump  discontinuities. 

Removable  singularities. 

xa-type  singularities. 

All  are  strictly  validated. 

Interval  Collection  Management  Component: 

Ordered  list. 

None  discarded. 

The  method  of  this  paper  does  not  belong  to  the  large  family  of  10c  or  so  algorithms 
considered  by  Rice  because  of  the  use  of  interval  computation  and  automatic  differentiation, 
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which  were  not  considered  in  i23>.  Details  of  the  actual  implementation  of  the  algorithms 
presented  here  in  an  environment  which  supports  interval  computation  will  be  given  in  §7. 
The  next  few  sections  describe  the  underlying  methodology. 


2.  Self-validating  Evaluation  of  Quadrature  Formulas 

Self-validation  of  numerical  computations  is  one  of  the  basic  motivations  of  interval  analysis 
[1],  [15],  [16].  The  goal  is  to  obtain  an  interval  which  contains  the  desired  result,  be  it 
real  or  set- valued.  In  the  case  of  the  integration  problem  (1.1),  a  self- validating  interval 
method  produces  an  interval  J  =  [c,  d\  which  is  guaranteed  to  contain  the  value  If  of  the 
integral.  The  width  of  this  interval  inclusion  will  depend  on  uncertainties  in  the  values  of 
the  integrand  and  the  limits  of  integration,  the  roundoff  error  in  the  actual  computation, 
and  the  truncation  error  appropriate  to  the  method  used.  All  of  these  quantities  can  be 
estimated  in  a  tedious  way  by  the  techniques  of  classical  error  analysis,  an  effort  which  is 
unnecessary  in  the  computational  environment  described  below.  However,  once  an  interval 
[e,d]  containing  I  is  found  by  whatever  method,  one  has  the  following  approximations  to 
I  and  corresponding  error  bounds  [21]: 

(21)  /*=j(«  +  <0,  i//-n< 

for  absolute  error,  or 


c  +  d' 


\If-l' 

U 


j  d  —  c 
[c+d 


for  relative  error,  with  cd  >  0  in  this  case.  It  follows  that  (1.2)  will  be  satisfied  if  an  interval 
J  =  [c,d]  can  be  obtained  with  width  w(J)  —  d  -  c  small  enough  so  that  w(J )  <  2a  and 
w(J)  <  P\c  +  d|. 

First,  the  problem  of  finding  an  interval  inclusion  J  of  //  will  be  considered.  The 
basic  method  for  interval  integration  by  use  of  standard  formulas  for  numerical  quadrature 
or  Taylor  series  was  first  described  by  Moore  i  1 4  .  To  illustrate  Moore’s  idea,  consider  a 
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standard  interpolator  integration  formula  of  the  form 

fb  A  fip)U)hr' 

(2.3)  /  f(x)dx  =  ^  wtf(xi)  -+-  cnh  - - - - , 

*'a  »=i  P‘ 

where  h  =  (6  -  a)/n,  and  a  <  £  <  b.  A  formula  of  type  (2.3)  will  be  called  a  quadrature 
formula  of  order  p  onn  points.  The  ordinary  Gauss  and  Newton-Cotes  integration  formulas 
follow  this  pattern  [7]. 

It  should  be  noted  that  integration  formulas  such  as  (2.3)  give  the  exact  value  of  the 
integral  If  of  functions  which  are  differentiable  p  times.  The  only  difficulty  is  that  the 
value  of  £  is  unknown.  For  practical  computation,  it  is  thus  customary  to  express  (2.3)  as 
the  sum  of  a  rule 

n 

(2.4)  rnf  =  ^2wtf(xi) 

i=i 

of  numerical  integration,  and  a  (truncation)  error  term 

(2.5)  enf  =  cnh-fp(t,h), 
where 

fP(t,h)  = 

denotes  the  Taylor  coefficient  of  order  p  in  the  expansion  of  /(£  +  h).  It  is  usual  to 
compute  1*  =  rnf  to  approximate  the  value  of  the  integral,  and  to  estimate  enf  somehow. 
Of  course,  if  /  is  a  polynomial  of  degree  p  -  1  or  less,  then  enf  =  0,  and  If  -  rnf. 

A  self-validating  computation  of  the  rule  rnf  of  numerical  integration  is  straightfor¬ 
ward  in  an  environment  which  provides  interval  arithmetic  and  monotone  interval  inclu¬ 
sions  of  the  library  functions  used  in  the  evaluation  of  f(x)  for  a  given  x.  Let  S  denote  the 
screen  of  floating-point  numbers  available,  and  IS  the  corresponding  set  of  closed  intervals 
ju,vj,  u,v  £  S.  If  /  is  evaluated  on  an  interval  A'  €  IS  using  interval  arithmetic  and 


fiF)mp 


\ 


library  functions,  then  the  result  is  the  natural  interval  inclusion  F(X)  of  /  on  X  such 
that 

(2.7)  f(X)={f(x)  ,T€  X}C  F[X). 

[15],  1 16].  If  W,,  X,  respectively  denote  the  smallest  intervals  in  IS  which  contain  the  real 
numbers  Wi ,  xt,  that  is,  IT,  —  [Vu/i,  Au>,-],  (Vzt,  Ax,],  where  V,A  denote  the  monotone 
downward  and  upward  roundings  from  the  real  numbers  R  to  S  { 1 1  j ,  then  the  inclusion 

n 

(2.8)  rnf  6  Rnf  =  '£wtF(Xt) 

i=  1 

is  guaranteed,  and  the  computation  of  Rn  can  be  done  automatically. 

An  automatic,  self-validating  computation  of  the  error  term  (2.5)  requires  an  addi¬ 
tional  ingredient.  This  consists  of  subroutines  for  the  generation  of  the  Taylor  coefficients 
/p(x,  h)  of  /.  These  use  well-known  recurrence  relations  for  the  arithmetic  operations  and 
library  functions  used  to  evaluate  f(x)  for  given  x  [6] ,  [15],  [16],  [19].  A  suitable  compu¬ 
tational  environment  provides  these  routines.  Corliss  and  Chang  |5]  have  shown  that  the 
calculation  of  /o(x,h),  /i(x,/i),...,/p(x,h)  requires  about 

(2.9)  t  =  ap2  +  bp+  c 

units  of  time,  where  a  depends  on  the  number  of  multiplications,  divisons,  and  calls  to 
library  functions  in  the  computation  of  /,  and  6  depends  on  the  number  of  additions  and 
subtractions  required.  In  any  case,  interval  evaluation  of  exactly  the  same  recurrence 
relations  yields  the  corresponding  interval  inclusions  Fq(X ,  H)  ,...,FP(X ,  H)  such  that 

(2.10)  fk{x,h)  <E  Fk{X,H),  k  =  0,1 . p, 

for  all  x  £  X  and  h  <E  H  [15],  [16],  Thus,  the  desired  interval  inclusion 

(2.11)  enfeEnf  =  CnH-Fr(X,H) 
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can  be  computed  automatically,  given  intervals  Cn.  H.X  6  IS  such  that  c„  6  Cn ,  h  6  //, 
and  ia.  6i  C  X.  (It  is  assumed  that  a  <  6;  the  contrary  case  can  be  handled  easily.)  It 
follows  from  the  above  that 

(2.12)  If  £  Hnf  ■+  Enf  ~  |c,dj, 

a  formula  for  automatic,  self-validating  numerical  quadrature.  Once  again,  if  /  is  a  poly¬ 
nomial  of  degree  p  -  1  or  less,  then  FP(X,H)  =  [0,0],  so  that  If  €  ZB„/,  and  the  width  of 
|c,d]  depends  only  on  the  roundoff  error  in  the  calculation  of  rnf. 

This  approach  was  the  basis  of  an  actual  computer  program  [9j,  which  met  the  accu¬ 
racy  criteria  (1.2),  if  possible,  by  choosing  H  sufficiently  small  [18].  However,  the  problem 
of  efficiency  was  not  addressed. 

Instead  of  splitting  the  numerical  integration  formula  (2.3)  into  a  rule  of  numerical 
integration  (2.4)  and  an  error  term  (2.5),  it  will  be  helpful  later  to  consider  it  to  represent 
the  scalar  product  of  the  augmented  function-value  vector 

(2-13)  f  =  (/(*i),..., /(**), /PU,fc))^ 

and  the  augmented  weight  vector 

(2.14)  w  =  (wu...,wn,  c„h), 

which  is  independent  of  the  integrand  /  and  depends  only  on  the  specific  formula  (2.3) 
used.  Thus, 

(2.15)  /  =  wf€W-F  =  J, 
where 

(2-16)  F  =  (F(X1),...,F(Xn),F;,(X,I/)) 

and 

(2-17)  W=  {Wu...,Wn.CnH) 

6 


are  the  corresponding  interval  inclusions  of  f,  w.  This  allows  the  computation  to  use  re¬ 
cently  developed  methods  for  highly  accurate  calculation  of  real  and  interval  scalar  prod¬ 
ucts  [12].  This  results  in  a  considerable  decrease  in  width  due  to  roundoff  error  in  the 
computed  value  of  J. 

The  integration  formula  (2.3)  can  be  interpreted  as  a  single-panel  rule,  or  as  a  multi¬ 
panel  rule,  meaning  that  a  simpler  formula  on  m  points  is  applied  k  times  to  the  corre¬ 
sponding  number  of  subintervals  of  X  —  [a,  6],  with  n  <  km.  Denoting  the  subintervals  of 
X  by  X{,  i  =  1,2 ,...,&,  this  means  that  an  integration  formula 


(2.18) 


f  m 

/  f{x)dx  =  Yl  Vijfixij)  + 
J Xi  j=l 


holds  in  each  subinterval,  where  hi  =  w(X{).  It  has  been  shown  [18]  that 

(2.19)  Vc^^X,)/^**)  •  C  cnw(X)F^(X)  • 

'  p!  p! 

In  addition  to  the  decrease  in  width  of  Fp{X,w(X))  by  a  factor  of  w(X)p  as  tu(X)  becomes 
small,  the  width  of  F(p)(X)  will  overestimate  the  width  of  /(p)(X)  by  less  as  u>(X)  — ►  0 
for  f^(x)  continuous  [15].  Thus,  the  gain  in  calculating  the  error  terms  over  smaller 
subintervals  can  be  substantial.  Roundoff  error  in  adding  a  number  of  interval  inclusions 
of  one-panel  rules  (2.18)  can  again  be  reduced  considerably  by  expressing  the  result  as  the 
scalar  product  of  the  extended  augmented  function-value  vector 

(2.20)  F  =  (F(Xn,...,F(Xlm),Fp(Xl,H1) . F(Xkl) . F(Xkm),  Fr(Xk,  Hk)) 


with  the  extended  augmented  weight  vector 


(2.21) 


W  =  (wn,  ...,wlm,ClmH1....,wk] . wkm,CkmHk). 


.  Taylor  Series  Methods 

he  seminal  paper  by  Moore  (14  j  also  provides  the  basis  for  self-validating  numerical  inte- 

\ 

ation  by  the  use  of  Taylor  series,  although  the  techniques  presented  by  him  in  this  case 
e  directed  toward  the  solution  of  the  initial-value  problem  for  ordinary  differential  equa- 
ans.  For  numerical  integration,  Taylor  series  are  more  appropriate  than  fixed  quadrature 
rmulas  for  interval-valued  endpoints  of  intervals  of  integration,  as  will  be  discussed  in  §6. 
irthermore,  Taylor  series  support  the  rigorous  approach  to  automatic  recognition  and 
eatment  of  singularities  described  in  §5. 

Of  course,  one  could  consider  (l.l)  to  be  the  solution  If  =  y(b)  of  the  initial-value 
oblem 


y'{x)  =  /(x),  y(a)  =  0, 


d  apply  Moore’s  methods  directly.  However,  since  f(x)  in  (3.1)  is  independent  of  y, 
ilike  the  usual  case  in  differential  equations,  it  is  simpler  to  use  the  capability  to  generate 
segment  of  the  Taylor  series  and  the  interval  remainder  term  automatically  to  perform 
ielf-validating  calculation  of  the  desired  integral. 

In  particular,  instead  of  expanding  the  solution  y(x)  of  (3.1)  at  x  =  a  as  in  the  case 
a  differential  equation,  it  is  advantageous  to  expand  f(x)  at  the  midpoint  c  =  [a  +  b)/2 
the  interval  X  =  [a,  6]  of  integration.  It  will  be  assumed  that  the  integrand  /  has  p  >  0 
rivatives  in  the  interval  of  integration.  For  h  =  (b  -  a)/2,  one  has 

\2  _  „\n-  1 

-+•••  +  /(n-1}(c) 


/(*)  =  /(c)  +  f'(c)[x  -c)  +  f'\c)  *  •  •  •  4.  ft-  (I  “  C) 


2! 


(»-!)! 


2) 


+  /<">(0^. 

n  <  p,  where  ^  €  X  is  between  c  and  x.  Let  F ^  be  an  interval  inclusion  of  on 
Then,  f^{£)  6  Fl'n^(X),  which  is  an  interval-valued  constant.  From  (3.2), 


1)  /(*)  6  f{c)  +  /'(c)(x  -  c) 


f{ 


(— i»(c)^ - 1 


(x-  c)(n- 

>  :  1)! 


F{n)(X) 
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Let 


(3.4) 


»(*)  =  /M(*  -  C)  +  /'(C)  +  •••  + 


n: 


be  an  indefinite  integral  of  the  Taylor  polynomial  of  degree  n  -  1  of  /(x).  Then, 

"t  4,  l-  -\n+  1 


(3.5) 

where 


f°  <t>  It  -  <«ln4 

/  /w^f»(x)  +e'"'(x)iT— i_ 

ya  la  ( n  +  lj! 


ir  ^ni 


n_1  i/t+1 

E  f",(c>(7TIj!  + 

« — n  v  1 


1=0 

t  even 


(3.6) 


+ 


rrn  +  1  ljn  +  l 

F(n)w#r TTT-^(n)W 


2F^n)(X) 


(»+!)! 

Hn41 


(»+!)! 


for  n  odd, 
for  n  even. 


(»  +  1)!’ 

Note  that  subtraction  does  not  “cancel”  equal  intervals  in  general.  One  has  [u,t>]  -  [u,  v]  = 

| u  -  v,v  -  u]  ^  [0,0]  unless  u  =  v,  in  which  case  the  interval  consists  of  a  single  point  [l], 
[15],  [16].  If  the  series  were  expanded  at  x  =  a  instead  of  x  =  c,  then  the  width  of  the 
interval  remainder  terms  would  be  increased  by  a  factor  of  2n+1. 

Formulas  (3.5)-(3.6)  resemble  (2.12)  for  ordinary  quadrature  rules,  with  the  evaluation 
of  the  integrand  at  n  points  replaced  by  its  value  and  the  values  of  its  first  n  -  1  derivatives 
at  a  single  point.  If  /  is  a  polynomial  of  degree  n  -  1  or  less,  then  F^(X)  =  [0,0],  and 
only  roundoff  error  effects  the  width  of  Jn. 

Jn  can  be  computed  directly  from  the  automatically  generated  interval  Taylor  coeffi¬ 
cients  Fk{C ,  H)  and  Fk(X,  H)  of  /,  k  =  0. 1 , . . . , n  <  p. 

Since 


(3.7) 


If 


[’  f(*)d 

J  a 


X  £  J  n  i 


n  =  0, 1 - ,p. 


the  intersection  principle  [8,,  [9j  can  be  used.  One  has 


(3.8) 


//e  n 


e  interval-valued.  One  strategy  for  handling  interval-valued  endpoints  is  to  allow  the 
icertainties  in  the  locations  of  the  endpoints  to  be  carried  over  into  uncertainties  in  the 
cations  of  the  nodes,  as  in  1NTE  ]9] .  For  example,  using  a  one-panel  Simpson’s  rule  on 
decimal  digit  machine  gives 

/'  ’  *  f{x)dx  €  — -r— ^  (F(j0,0.1j  +  4F(|l.55,1.65])  +  F(|3. 1,3.2)) 

J  [0,0.1]  b 

1.7) 

here  A  =  [0,0.1],  B  =  [3.1, 3.2).  For  f(x)  =  sini,  (6.7)  yields  [1.880,2.214]  using  the 
ccurate  scalar  product,  while  the  correct  answer  is  [1.994,2.000]. 

A  better  strategy  is  to  concentrate  the  uncertainty  at  the  ends: 

y  (3.11,(3.21  rO.l  [3.1  /-  [3. 1 ,3.2] 

3.8)  /  f(x)dx=  I  f(x)dx+  /  f(x)dx+  I  f(x)dx. 

J\ 0,0.1]  •'10,0.1]  J  0.1  j  3.1 

’he  middle  part  is  handled  as  an  RR  integral  as  described  in  §§3-4,  while  the  other  two 
itegrals,  of  the  types  IR  and  RI,  respectively,  will  be  treated  in  the  manner  to  be  described 
ielow.  In  this  example,  we  can  get  j  1.984 ,2.073].  As  the  widths  of  the  endpoints  increase, 
he  advantage  of  the  second  strategy  becomes  more  pronounced. 

The  general  case  (6.3)  can  be  expressed  in  terms  of  integrals  of  the  above  types.  This 
ase  in  turn  splits  into  six  subcases,  according  as  A  and  B  are  (i)  disjoint,  (ii)  overlapping, 
•r  (iii)  one  is  contained  in  the  other.  More  precisely,  for  A  =  [  A  j, ,  .4  ft  j ,  B  =  [Bl,Br\, 
uppose  that  Ar  <  Br.  Then  we  have: 


Case  1.  AR  <  BL  (disjoint  interval  endpoints).  Here,  the  integral  (6.3)  can  be  written 


fti  r An  rtlL  r\BL,Br] 

6.9)  /  f(x)dx=  /  f(x)dx  +  /  f(x)dx  +  /  f(x)dx , 

Ja  J\Al,Ap\  j  Ap  J  Bl 

ind  thus  is  the  sum  of  integrals  of  types  IR,  RR,  and  Rl.  respectively. 

Case  2.  AL  <  BL  <  AR  <  BR  (overlapping  interval  endpoints).  Here, 


r  B  r  B  l  f  Ri'Ai.  r\Ap,tiL, 

/  f{x)dx  =  /  f(x)dx  -  /  f[r)dx  -  /  f{x)dx, 

J  A  J  \Ai  .Bl  i  J  B I  .Ar.  *  A  p 


\Ap  ,Bl 


•V-'.v.-.-.v.  . 


I*.  c_.-_ 


>  .-. . 
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ie  natural  interval  inclusions  Fk{X.C] _ ,C„, )  of  /  and  its  Taylor  coefficients  on  an 

erval  A'  €  IS  are  again  obtainable  on  a  computer  by  using  interval  computation  and 
tomatic  differentiation.  In  particular,  for  an  interval  polynomial  (6.2),  FP(X,H)  =  [0, 0] 
p  >  m,  just  as  in  the  real  case.  The  definition 

4)  /(x,Ci,...,  Cm)dx  —  |  J  -/"(j,ci,...,cfri)dz  j  c  j  G  C\ , . . . ,  cm  G 

cribes  the  type  of  integrals  to  which  the  methods  of  this  paper  apply.  It  is  assumed  that 
B  and  the  coefficient  intervals  Ct  all  belong  to  IS,  and  hence  are  machine-representable, 
lecessary,  outward  rounding  can  be  used  to  obtain  them  from  real  intervals,  preserving 
lusion  of  the  desired  integral. 

On  the  basis  of  (2.3),  it  is  to  be  expected  that  the  best  results  will  be  obtained  for 
!grands  /  which  are  very  smooth  as  functions  of  x.  However,  one  can  always  compute 
interval  Riemann  sum  F(Ar)  •  w[X).  This  is  self-validating,  because 

.)  f  f{x)dx  C  F(X)  ■  w(X), 

Jx 

is  inaccurate  and  slow  to  converge  [4],  [20].  It  is  important  to  be  able  to  confine  bad 
avior  of  the  integrand  to  very  small  intervals  for  (6.5)  to  be  useful. 

The  meaning  of  tolerance  for  problems  with  interval-valued  endpoints  requires  some 
ification.  Consider  the  problem 


i  compute  J  =  [-0.004,0.79],  for  example,  then  w(J )  =  0.794.  although  the  estimate 
ach  endpoint  is  in  error  by  less  than  0.005.  Hence,  a  requested  tolerance  must  be  large 
gh  to  accomodate  the  uncertainties  which  are  inherent  in  the  problem  being  solved. 
The  cases  that  one  or  both  endpoints  of  integration  are  nondegenerate  intervals  in 
dll  now  be  considered.  The  possible  situations  will  be  denoted  by  IR.  Rl.  and  II. 
‘Ctively.  according  as  the  lower,  upper,  or  both  endpoints  of  the  interval  of  integration 
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6.  Extension  to  Interval  Values 


The  definition  of  the  integral  (1.1 )  has  been  extended  to  arbitrary  interval-valued  inte¬ 
grands  /  [4],  1 20].  For  smooth,  real-valued  functions  such  as  those  considered  in  §2,  the 
interval  integral  and  the  Riemann  integral  coincide.  The  concept  of  interval  integration 
is  useful  in  connection  with  integrands  which  have  jump  discontinuities  or  singularities  of 
various  kinds.  The  definition  given  in  |4|  can  be  extended  to  the  case  that  the  limits  of 
integration  are  interval-valued: 

for  A,B  €  IR,  the  set  of  all  finite  intervals  with  real  endpoints. 

There  are  several  reasons  to  want  to  be  able  handle  interval  endpoints  of  integration. 
First  of  all,  some  real  numbers,  such  as  n,  cannot  be  represented  exactly  in  S,  and  have 
to  be  replaced  by  the  corresponding  small  intervals  such  as  [Vtt,  Act]  in  IS.  Secondly,  the 
limits  of  integration  may  come  from  measurements  other  estimates,  and  thus  are  known 
to  lie  between  certain  limits  even  though  their  exact  values  are  uncertain.  Finally,  one  can 
be  interested  in  bounds  for  the  value  of  the  integral  (l.l)  over  ranges  of  values  of  limits 
of  integration.  In  this  case,  rather  than  satisfy  an  accuracy  criterion  such  as  (1.2),  it  is 
usually  desired  to  find  an  interval  J  containing  (6.1)  which  is  as  small  as  possible. 

Similar  considerations  apply  to  interval-valued  integrands.  The  case  which  usually 
arises  in  practice  is  that  the  integrand  /  is  a  function  not  only  of  the  independent  variable 
x,  but  also  several  parameters  cj,C2,. . .  ,cm.  For  example,  /  could  be  a  polynomial  of 
degree  m  -  1  with  coefficients  determined  by  observations, 

(6.2)  /(x)  =  C,  +  C2x+---  +  Cmxm-1,  C{  €  IS,  i  —  1,2, . . . ,  m. 

In  general,  given  intervals  C\,  C2, . . .  ,  Cm,  it  is  natural  to  define 

/(x;  C\. . . .  ,Cm)  =  {/(x; c j , . . . , cm)  I  Cj  €  C\ ,...,cm  £  Cm}. 


f{x)dx 


£  A,  be  B 


}• 


(6.1) 


/ 


B 


f(x)dx 


(6.3) 


then  the  integrand  cannot  be  handled  by  this  method,  and  a  message  to  that  effect  is  sent 
to  the  user. 

These  methods  for  recognizing  and  handling  certain  singularities  extend  the  domain 
of  applicability  of  the  program  to  include  many  integrands  which  arise  in  applications. 
However,  their  usefulness  is  somewhat  limited.  For  one  thing,  the  location  of  the  singu¬ 
larity  must  be  known  in  advance,  or  guessed.  For  popular  sets  of  test  problems,  guessing 
one  or  both  endpoints  usually  works.  For  real  problems,  locations  of  singularities  may  be 
unknown.  However,  the  method  given  here  validates  correct  guesses.  The  endpoints  to  be 
investigated  for  the  presence  of  singularities  must  be  machine  numbers,  not  intervals.  If 
an  interval-valued  endpoint  is  even  one  machine  number  wide,  then  the  problem  contains 
integrands  which  are  unbounded  on  a  set  of  positive  measure.  If  a  singularity  is  in  the 
interior  of  the  interval  of  integration,  then  we  can  determine  its  location  to  within  one  or 
two  machine  numbers.  From  this,  its  contribution  to  If  could  be  estimated,  but  the  possi¬ 
bility  that  the  integrand  is  unbounded  on  a  set  of  positive  measure  cannot  be  eliminated. 
Since  a  validated  answer  cannot  be  produced  in  this  case,  an  error  return  is  selected.  If 
the  user  knows  that  the  singularity  occurs  at  a  machine  number,  then  the  integral  should 
be  calculated  as  the  sum  of  two  integrals,  with  the  singularity  at  one  endpoint  of  each. 


Suppose  that  the  integrand  has  the  form 


(5.9) 


f{x)  =  (a  -  x)  *<t>(x), 


where  4>(x)  is  analytic  at  x  =  a.  If  c  is  chosen  closer  to  a  than  any  other  singularity,  then 
the  series  for  /  expanded  at  x  =  a  is  asymptotic  to  the  series  for  v(i)  =  (a  -  x)“*.  The 
Taylor  coefficients  u,  =  ut(a,  h)  of  v  satisfy  the  recurrence  relation 


(5.10) 


v.+i  =  Vi 


s  -  1\  h 

)WC' 


where  Rc  is  the  radius  of  convergence  of  the  series.  If  /  cannot  be  evaluated  on  (a,  6],  then 
one  attempts  to  find  constants  Ki ,  Hr,  $l,  sr  such  that 


(5.11)  KL(a  -  x)~“  <  f(x)  <  Kp(a  - 

for  x  near  a.  If  such  constants  can  be  found  and  (5.11)  validated,  then  we  proceed.  If 
si  >  1,  then  the  program  can  guarantee  that  If  does  not  exist.  We  know  of  no  other 
numerical  quadrature  routine  which  can  validate  nonexistence  of  an  integral.  If  sr  <  1, 
then 


(5.12)  -^-(a-a')1-^  <  f  f(x)dx  <  -^-(c  -  a')1^- 

1  _  Ja  1  -  SR 

These  bounds  can  be  made  as  tight  as  desired  by  taking  a'  close  enough  to  a.  The  interval 
a,  a' j  is  placed  on  the  list  of  subintervals,  and  processing  continues  on  the  subinterval 
a\ br 

A  singularity  at  b  can  be  handled  similarly. 

If  Kl,  Kr ,  si,  sr  cannot  be  found,  or  if 

(5.13)  sL  <  1  <  Sr , 
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because  of  its  width.  All  other  subintervals  can  be  processed  using  higher-order  rules.  This 
strategy  applies  to  singularities  in  f  which  occur  anywhere  in  the  interval  of  integration, 
and  not  just  at  endpoints.  It  also  works  for  singularities  in  higher  derivatives,  whose 
presence  might  not  even  be  known  to  the  user.  Once  a  singularity  of  this  type  has  been 
confined  to  a  sufficiently  small  subinterval,  it  is  possible  to  meet  reasonable  requirements 
for  accuracy,  with  self- validation. 

5.2.  Jump  discontinuities 

If  the  integrand  is  given  by  an  ordinary  mathematical  expression,  then  it  is  difficult  to 
represent  a  function  which  has  jump  discontinuities,  and  yet  can  be  evaluated  at  every  point 
in  the  interval  of  integration.  However,  the  user  can  supply  a  subroutine  for  evaluation 
of  the  integrand  which  produces  jump  discontinuities.  These  appear  to  the  program  as 
singularities  in  /'.  On  any  subinterval  which  contains  a  jump,  only  the  Riemann  sum 
is  available,  and  its  width  will  lead  to  frequent  selection  of  this  subinterval  for  further 
processing,  as  before.  No  special  algorithms  are  needed  in  this  case. 

This  behavior  is  similar  to  CADRE  [b2j.  Upon  recognizing  a  jump,  CADRE  subdi¬ 
vides  the  interval  and  uses  a  low-order  rule. 

5.3.  Some  removable  singularities 

The  Taylor  series  method  permits  handling  of  some  removable  singularities.  Consider 
the  problem 

(5-8)  //=  -dx, 

Jo  x 

in  which  the  integrand  has  a  removable  singularity  at  x  —  0.  In  this  case,  the  integrand 
cannot  be  evaluated  directly  on  any  interval  containing  0,  but  /  can  be  expanded  at  x  =  0 
using  I’Hospital’s  rule,  which  can  be  applied  automatically  [6j.  A  short  Taylor  series  with 
remainder  for  /  at  0  is  sufficient  to  bound  If  near  0,  and  the  rest  of  the  interval  of 
integration  is  processed  in  the  normal  manner. 

5.4.  Some  algebraic  singularities  at  endpoints 
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Although  the  original  problem  (5.2)  is  well-behaved,  the  interval  integral  in  (5.3)  contains 
integrals  such  as 


(5.4) 


-  x  dx , 


and  consequently  does  not  exist.  Hence,  the  correct  response  is  that  (5.3)  cannot  be 
evaluated  on  the  entire  interval  of  integration,  just  on  jO,  Vx],  for  example.  This  suggests 
that  standard  numerical  quadrature  routines,  which  avoid  evaluation  of  the  integrand  at 
endpoints,  can  be  fooled  into  returning  values  for  integrals  such  as 


for  c  sufficiently  small,  instead  of  informing  the  user  of  possible  difficulty.  The  types  of 
singularities  which  can  be  handled  by  the  techniques  presented  here  will  now  be  described. 

5.1.  Singularities  in  derivatives  of  /. 


The  integrand  /  may  not  have  enough  derivatives  for  some  rules  which  the  integration 
program  can  apply.  In  the  process  of  automatic  generation  of  interval  Taylor  coefficients 
on  an  interval  X,  the  nonexistence  of  a  derivatives  of  /  beyond  a  certain  order  is  detected 
and  reported.  As  long  as  /  itself  can  be  evaluated  on  the  entire  interval  of  integration, 
the  order  adaptation  strategy  handles  singularities  in  the  derivatives  of  /.  For  example, 
consider  the  problem 


(5.6) 


//=  /  y/idx=  1 
Jo 


The  first  derivative  f  is  undefined  at  x  =  0,  so  the  only  rule  that  can  be  applied  to  the 
entire  interval  of  integration  is  the  Riemann  sum 


(5.7) 


I  y/xdx  €  min  v/x,  max  v/x  =  (0, 8] 

Jo  Mo, 4)  10,4)  J  1  J 


At  this  point  the  subinterval  adaptive  strategy  takes  over.  At  each  step,  the  subinterval 
containing  0  requires  a  Riemann  sum,  and  thus  is  frequently  selected  for  further  processing 


5.  Treatment  of  Singularities 

Among  the  quadrature  routines  which  provide  estimates  for  If,  the  more  successful  ones 
have  special  provisions  to  handle  and  perhaps  recognize  certain  types  of  singularities. 
QUADPACK  [17],  for  example,  can  handle  integrals  of  the  form 

I1, 

(5.1)  I  (x  -  a)a(b  -  x)^v(x)f(x)dx, 

J  a 

where  a, (3  >  -1  and  v(x)  =  1,  log(x  -  a),  log(6  -  x),  or  log(x  -  a)  log(6  -  x),  provided 
the  user  supplies  a ,  /?,  and  the  form  of  v.  CADRE  [3]  attempts  to  detect  and  verify  the 
presence  of  jump  discontinuities  or  x^-type  singularities  in  the  integrand.  The  program 
described  here  attempts  to  recognize  and  handle  automatically 

(1)  singularities  in  derivatives  of  /, 

(2)  some  removable  singularities, 

(3)  jump  discontinuities,  and 

(4)  some  algebraic  singularities  at  endpoints. 

Since  guaranteed  bounds  are  computed  for  If,  the  task  of  such  a  program  is  more 
difficult  than  for  methods  which  yield  only  an  estimate.  For  example,  Gaussian  quadrature 
can  be  used  in  the  neighborhood  of  an  endpoint  singularity,  because  the  integrand  need 
not  be  evaluated  at  the  endpoint.  To  give  guaranteed  bounds,  however,  requires  that  the 
function,  or  some  of  its  derivatives,  be  evaluated  on  the  entire  interval.  Hence,  the  price 
of  the  guarantee  is  a  restriction  on  the  applicability  of  the  program.  Either  it  verifies  that 
the  integrand  is  in  its  domain  of  applicability,  or,  if  it  cannot  deliver  guaranteed  bounds, 
it  notifies  the  user  with  an  indication  of  the  difficulty. 

An  example  will  indicate  what  can  go  wrong  at  endpoints.  For  example, 


could  be  presented  to  the  computer  as 

/•  |  Vir,  An-j 

(5.3)  I2f  -  I  v/l^77r5  Ajt]  -  x  dx. 

Jo 
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Because  of  the  difference  between  the  remainder  terms  in  (3.6)  for  odd  and  even  n,  it 
is  prudent  to  calculate  at  least  one  extra  value  of  J„.  However,  the  conditions  above 
guarantee  that  no  more  than  two  series  terms  will  be  computed  beyond  the  one  which 
gives  the  narrowest 

Another  important  strategy  for  order  adaptation  involves  the  case  that  the  integrand 
has  singularities  in  a  certain  derivative  in  the  given  subinterval.  This  will  be  discussed 
more  fully  in  the  next  section.  If  a  certain  derivative  cannot  be  evaluated,  this  will  be 
detected,  and  the  method  will  be  restricted  to  rules  or  orders  of  Taylor  expansion  which 
use  only  the  derivatives  which  can  be  evaluated. 

The  strategy  for  subinterval  adaptation  retains  ail  subintervals.  At  each  step,  the 
subinterval  which  makes  the  largest  contribution  to  the  width  of  J  is  processed  by  breaking 
it  into  further  subintervals.  This  processing  continues  until 

(i)  vu(J)  is  small  enough  to  satisfy  the  accuracy  requirement  (1.2), 

(ii)  the  noise  inherent  in  function  evaluation  limits  further  reduction  of  w(J),  or 

(iii)  more  than  the  maximum  number  M  of  function  evaluations  have  been  performed. 

The  second  termination  criterion  in  this  list  is  particularly  important.  If  the  noise  in  the 
function  evaluation  is  large  relative  to  the  accuracy  requested,  eventually  the  width  of  the 
truncation  error  Ef  is  made  so  small  that  it  adds  nothing  to  the  width  of  the  rule  Rf 
alone.  At  this  point,  further  increase  in  accuracy  is  not  possible.  Without  the  guaranteed 
bounds  provided  by  the  interval  computation,  many  standard  methods  cannot  recognize 
when  this  point  has  been  reached,  and  the  calculation  should  be  terminated.  Malcolm 
and  Simpson  j  13]  observe  that  the  strategy  of  processing  the  worst  subinterval  results 
in  local  errors  of  roughly  equal  magnitude.  Rice  [23]  calls  this  an  ordered  list  interval 
collection  management  component,  and  lists  some  advantages  and  disadvantages  which 
will  be  discussed  in  more  detail  in  §7. 


Alternatively,  the  actual  approximate  integrals  J,  could  be  examined,  but  this  requires 
more  computation  than  use  of  the  error  terms  alone.  Suppose  that  p  denotes  the  maximum 
value  of  the  width  u;(F(C))  of  the  interval  evaluation  of  /  at  a  node  C  of  the  integration 
formula  being  used.  The  total  cost  can  be  taken  to  be  proportional  to  n,  •  u>( J,),  where  n, 
is  the  number  of  function  evaluations.  The  width  w(Jt)  can  be  estimated  by  p  -+  w{EnJ)> 
and  thus  i  can  be  chosen  as  the  minimizer  of 

(4.2)  m(t)  =  min{ty(Jo),ni(M  +  #«,/)}■ 

Thus,  more  nodes  are  used  in  a  given  subinterval  only  if  a  significant  reduction  in  width 
of  the  approximate  integral  results.  It  should  also  be  noted  that  a  given  integration  rule 
can  be  used  in  several  integration  formulas  having  remainder  terms  of  different  orders. 
For  example,  Stroud  and  Secrest  [24]  give  error  terms  for  Gaussian  integration  rules  on  n 
nodes  which  have  orders  1  <  p  <  2 n.  In  certain  cases,  the  error  terms  corresponding  to 
smaller  than  maximum  p  can  be  narrower  (for  exam^u;  for  highly  oscillatory  integrands), 
and  the  use  of  these  formulas  instead  of  the  standard  ones  will  minimize  vu(Jn). 

In  the  case  of  Taylor  series,  the  intersection  (3.9)  procedure  can  result  in  approxima¬ 
tions  In  which  are  considerably  better  than  Jn  given  by  (3.6),  which  can  be  viewed  as  the 
sum  of  a  rule  and  an  error  term.  The  intervals  In  are  monotone  decreasing  in  width,  so 
the  increase  in  accuracy  has  to  be  balanced  against  the  cost  in  time  (2.9)  of  generating 
more  terms  of  the  series.  The  constants  a,b,  c  in  (2.9)  can  always  be  determined  for  a 
given  integrand,  so  the  corresponding  heuristic  can  be  based  on  the  function 

(4.3)  0(n)  =  u>(/„)  •  (a**2  +  bn  +  c). 

Generation  of  the  Taylor  series  can  be  stopped  when 

(i)  f n -  2  “  f n -  1  =  Ini 

(ii)  0(n  -  2)  <  0(n),  or 
(iii)  n  =  p. 
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4.  Adaptive  Strategies 

The  computer  program  INTE  (9j  demonstrated  the  reality  of  automatic,  self-validating 
numerical  quadrature  using  Newton-Cotes  and  Gaussian  integration  formulas  in  the  way 
discussed  in  §2.  (Euler-Maclaurin  integration  was  added  to  the  capabilities  of  INTE  later 
[10].)  However,  there  was  no  attempt  to  address  the  problem  of  efficiency,  a  defect  reme¬ 
died  in  the  program  described  later.  In  order  to  satisfy  the  accuracy  criterion  (1.2),  if 
possible,  INTE  simply  divided  the  interval  of  integration  into  a  sufficient  number  of  equal 
subintervals  [9],  [18].  By  contrast,  popular  numerical  integration  packages  such  as  CADRE 
[3]  and  QUADPACK  [17]  use  information  obtained  about  the  behavior  of  the  integrand  to 
attempt  to  reduce  the  number  of  function  evaluations  to  a  minimum.  The  method  given 
here  uses  similar  strategies,  except  that  estimates  of  the  error  based  on  evaluation  of  the 
integrand  at  a  finite  set  of  points  are  replaced  by  guaranteed  bounds.  This  eliminates  the 
need  for  “safety  factors”  [23]. 

Adaptive  strategies  fall  into  the  categories  of  order  adaptation,  which  relates  to  the 
choice  of  the  formula  used  in  each  subinterval,  and  subinterval  adaptation,  which  de¬ 
termines  how  the  original  interval  of  integration  is  broken  up  into  subintervals.  Order 
adaptation  is  somewhat  simpler  than  subinterval  adaptation,  and  will  be  considered  first. 
For  methods  based  either  on  standard  quadrature  formulas  or  Taylor  series,  order  zero 
refers  to  the  interval  Riemann  sum  F(X)  •  w(X),  which  always  contains  fx  f(x)dx  [4]. 

Given  a  suite  of  numerical  integration  formulas  of  the  form  (2.3),  suppose  that  X  €  IS 
is  the  current  subinterval  of  integration.  Specifically,  suppose  that  the  given  rules  are  of 
order  i  for  i  =  0, 1, . . . , k  on  »i  points.  Once  the  interval  Taylor  coefficients  F,(X ,  H)  have 
been  formed,  then  the  order  p  <  k  of  the  most  accurate  rule  can  be  chosen  to  be  the  value 
of  i  for  which  the  width  of  the  error  term 

(41)  EnJ  =  CniH  •  Fi(X,H) 

is  minimum,  i  =  0,1,..., k. 
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Alternatively,  the  actual  approximate  integrals  J,  could  be  examined,  but  this  requires 
more  computation  than  use  of  the  error  terms  alone.  Suppose  that  p  denotes  the  maximum 
value  of  the  width  w(F(C))  of  the  interval  evaluation  of  /  at  a  node  C  of  the  integration 
formula  being  used.  The  total  cost  can  be  taken  to  be  proportional  to  n,  •  tu(  Jt),  where  n, 
is  the  number  of  function  evaluations.  The  width  w(Jt)  can  be  estimated  by  p  +  w(EnJ ), 
and  thus  i  can  be  chosen  as  the  minimizer  of 


m(i)  =  min{te(</0),ni(/i  +  w(£n,/)}- 


Thus,  more  nodes  are  used  in  a  given  subinterval  only  if  a  significant  reduction  in  width 
of  the  approximate  integral  results.  It  should  also  be  noted  that  a  given  integration  rule 
can  be  used  in  several  integration  formulas  having  remainder  terms  of  different  orders. 
For  example,  Stroud  and  Secrest  [24]  give  error  terms  for  Gaussian  integration  rules  on  n 
nodes  which  have  orders  1  <  p  <  2 n.  In  certain  cases,  the  error  terms  corresponding  to 
smaller  than  maximum  p  can  be  narrower  (for  example,  for  highly  oscillatory  integrands), 
and  the  use  of  these  formulas  instead  of  the  standard  ones  will  minimize  w(Jn). 

In  the  case  of  Taylor  series,  the  intersection  (3.9)  procedure  can  result  in  approxima¬ 
tions  /„  which  are  considerably  better  than  Jn  given  by  (3.6),  which  can  be  viewed  as  the 
sum  of  a  rule  and  an  error  term.  The  intervals  /„  are  monotone  decreasing  in  width,  so 
the  increase  in  accuracy  has  to  be  balanced  against  the  cost  in  time  (2.9)  of  generating 
more  terms  of  the  series.  The  constants  a,b,c  in  (2.9)  can  always  be  determined  for  a 
given  integrand,  so  the  corresponding  heuristic  can  be  based  on  the  function 

(4.3)  0(n)  =  w(ln)  ■  (an2  +  bn  +  c). 

Generation  of  the  Taylor  series  can  be  stopped  when 
(*)  In  —  2  =  In-1  =  Ini 

(ii)  0(n  -  2)  <  6(n),  or 

(iii)  n  =  p. 


/•V*  .*■  «' 

.--V-V  V-  . 
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4.  Adaptive  Strategies 

The  computer  program  INTE  [9]  demonstrated  the  reality  of  automatic,  self-validating 
numerical  quadrature  using  Newton-Cotes  and  Gaussian  integration  formulas  in  the  way 
discussed  in  §2.  (Euler-Maclaurin  integration  was  added  to  the  capabilities  of  INTE  later 
(lOj.)  However,  there  was  no  attempt  to  address  the  problem  of  efficiency,  a  defect  reme¬ 
died  in  the  program  described  later.  In  order  to  satisfy  th  *»xcuracy  criterion  (1.2),  if 
possible,  INTE  simply  divided  the  interval  of  integration  into  a  sufficient  number  of  equal 
subintervals  [9],  [18).  By  contrast,  popular  numerical  integration  packages  such  as  CADRE 
[3]  and  QUADPACK  jl7j  use  information  obtained  about  the  behavior  of  the  integrand  to 
attempt  to  reduce  the  number  of  function  evaluations  to  a  minimum.  The  method  given 
here  uses  similar  strategies,  except  that  estimates  of  the  error  based  on  evaluation  of  the 
integrand  at  a  finite  set  of  points  are  replaced  by  guaranteed  bounds.  This  eliminates  the 
need  for  “safety  factors”  [23]. 

Adaptive  strategies  fall  into  the  categories  of  order  adaptation,  which  relates  to  the 
choice  of  the  formula  used  in  each  subinterval,  and  subinterval  adaptation,  which  de¬ 
termines  how  the  original  interval  of  integration  is  broken  up  into  subintervals.  Order 
adaptation  is  somewhat  simpler  than  subinterval  adaptation,  and  will  be  considered  first. 
For  methods  based  either  on  standard  quadrature  formulas  or  Taylor  series,  order  zero 
refers  to  the  interval  Riemann  sum  F(X)  ■  u>(AT),  which  always  contains  fx  f(z)dx  [4]. 

Given  a  suite  of  numerical  integration  formulas  of  the  form  (2.3),  suppose  that  X  G  IS 
is  the  current  subinterval  of  integration.  Specifically,  suppose  that  the  given  rules  are  of 
order  i  for  i  =  0, 1, . . . , k  on  n,  points.  Once  the  interval  Taylor  coefficients  F,(X,  H )  have 
been  formed,  then  the  order  p  <  k  of  the  most  accurate  rule  can  be  chosen  to  be  the  value 
of  t  for  which  the  width  of  the  error  term 

(4.1)  EnJ  =  Cn,H  Fi(X,H) 


is  minimum,  i  =  0, 1, . . . ,  k. 


This  intersection  can  be  calculated  as  the  corresponding  interval  Taylor  coefficients  of  the 
integrand  are  generated: 


(3.9) 


10  =  J0,  (a  Riemann  sum), 

In  ~  In  —  1  0  fl  —  1,2 , ,p. 


This  provides  a  means  to  determine  the  highest  useful  term  of  the  Taylor  expansion,  since 
the  calculation  can  be  terminated  when  effective  decrease  in  the  widths  of  the  intervals 
{/„}  ceases,  or  when  the  desired  tolerance  is  met. 

As  in  the  case  of  (2.12),  formula  (3.6)  can  be  evaluated  as  the  scalar  product  of  the 
vectors 


(3.10)  F  =  {Fo(C,H),F1{C,H),...,Fn-i{C,H),Tn{C,H)), 

where  Tn{C,H)  =  Fn{C,H)  -  Fn{C,H)  or  Tn(C,H)  =  2 Fn{C,H)  according  as  n  is  odd 
or  even,  and  the  vector 

(3.11)  W  =  (W0,Wu...,Wn-UWn), 
where 

(3.12)  Wk  =  k  =  0,1,. ..,n. 

The  width  of  Jn  can  also  be  reduced  somewhat  if  the  expansion  is  about  C  =  [c,cj, 
c  6  S.  In  this  case,  the  interval  stepsize  H  might  have  to  be  increased  a  very  small 
amount  to  maintain  inclusion  of  the  subinterval  of  integration.  In  §5,  it  will  be  noted  that 
expansion  about  an  endpoint,  as  in  (3.1),  or  some  other  point  in  the  interval  of  integration, 
can  be  helpful  if  the  integrand  has  removable  singularities. 
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the  sum  of  integrals  of  types  IR,  II,  and  RI. 


Case  3.  BL  <  AL  <  AR  <  BR  (A  is  properly  contained  in  B).  Here, 

r  B  rAi  r\At.,An]  r\Ap ,  Bp\ 

(6.11)  /  f(x)dx=  -  /  f(x)dx  -  /  f(x)dx+  /  f{x)dx, 

Ja  J\Bl,Al\  J \Al  ,Ap  ]  JAp 

again  the  sum  of  integrals  of  types  IR,  II,  and  RI. 


The  other  three  cases  are  obtained  for  Bn  <  Ap  by  reversing  the  roles  of  A  and  B  in  the 
preceding. 

Case  1  is  undoubtedly  the  one  which  is  encountered  most  often  in  practice.  To  illus¬ 
trate  Case  3,  consider 


r  [0,5]  f2  /“  1 2 ,3]  r\3,5] 

(6.12)  /  f(x)dx  =  ~  /  f(x)dx+  I  f(x)dx+  I  f(z)dx. 

7,2,3!  |0,2]  7,2, 3j  73 

We  are  not  aware  of  physical  problems  which  give  rise  to  integration  problems  other  than 
Case  1  (except  to  bound  the  values  of  integrals  over  ranges  of  endpoints),  but  provision 
for  these  cases  adds  little  to  the  machinery  which  is  required  to  handle  Case  1. 

Integrals  of  types  RI,  IR,  and  II  can  be  handled  by  Gauss  or  Newton-Cotes  formulas 
with  interval-valued  nodes,  but  this  leads  to  wide  interval  bounds.  The  method  of  Taylor 
series,  as  outlined  in  §4,  applies  as  easily  to  the  cases  of  one  or  both  endpoints  interval¬ 
valued  as  to  RR  type  integration.  Hence,  the  program  uses  Taylor  polynomials  for  integrals 
of  types  RI,  IR,  and  II,  even  if  the  user  has  chosen  Gauss  or  Newton-Cotes  formulas  for 
use  on  type  RR  subintervals. 

Using  the  same  notation  as  in  §3,  ^(z)  denotes  an  indefinite  integral  (3.4)  of  f(x). 
The  one-panel  form  of  the  various  integration  formulas  are  then: 

Type  RI: 


(6.13) 


L 


Ml 


f(x)dx  C  g{x) 


I0.*) 


F(n)(|a,  6]) 


(z  -  c)n+1  1MI 

(rc+*)!  |a 


Type  IR: 
(6.14) 


f(x)dx  C  g(x)  \  _  +7,")(;a.6])5- — 


|f. 

i 


Type  II: 


(6J5) 


+  F^(\a,b)) 


{x~c)n+1 
(n  +  1)! 


!«.*»] 
(a, 6] 


Each  uses  intersections  of  subsequent  estimates  and  order  adaptation  as  the  RR  algorithm 
does. 

Formulation  of  the  multipanel  forms  for  the  above  types  requires  some  care.  Suppose 
the  nodes  a  =  xq  <  xj  <  •••  <  x*  =  b  are  selected,  and  set  Yt  =  [*,•_!,*<]  for  i  = 
1,2,..., A:.  For  type  II  integrals,  let 


(6.16) 


The  subintervals  Yt  are  chosen  by  the  same  subinterval  adaptive  strategy  as  used  for  type 
RR  integrals.  Thus,  the  algorithm  for  type  II  integrals  is  the  same  as  for  type  RR.  except 
that  the  integral  on  each  subinterval  is  computed  using  the  one-panel  type  II  rule,  equation 
(6.15),  with  intersection. 
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The  multipanel  rules  for  types  RI  and  IR  do  not  lend  themselves  to  subinterval  adap¬ 


tation.  By  definition, 


(6.20) 


J  f{x)dx=^yj  f[x)dx  t  e  [a,  6]  J 

=  \j{f'f(x)dx  \teY'} 

fYl 

=  /  f(x)dx+ 

J  a 

fii  rY 2 

+  /  f(x)dx  +  /  f(x)dx+ 

j  a  J  X\ 


+  •••  + 


f  f(x)dx-\-  f  f(x)dx+---+  j  f(x)da 

J  a  J  X  i  1 


Once  the  endpoints  of  a  subinterval  are  chosen,  that  subinterval  can  be  subdivided  only 
with  great  difficulty.  We  observe  that 


(6.21) 


/■(a, 6]  r\a,  6] 

/  f{x)dx  c  j  f(x)dx, 
Ja  •/  }<x,6J 


so  the  type  II  algorithm  gives  bounds  using  adaptation.  If  those  bounds  are  not  small 
enough,  then  we  apply  the  algorithm  given  by  equation  (6.20)  using  a  fixed  stepsize  equal 
to  the  smallest  stepsize  chosen  adaptively  by  the  type  II  algorithm. 

Type  IR  integrals  are  done  similarly. 

In  the  general  case  of  the  integral  (6.3),  the  allowable  tolerances  and  the  maximum 
number  of  function  evaluations  must  be  apportioned  among  three  distinct,  independent 
integrations.  These  are  given  in  proportion  to  the  width  of  the  respective  subintervals.  If 
r(Xt)  denotes  the  tolerance  allowed  on  the  subinterval  X,,  r*  is  the  tolerance  remaining, 


(6.21) 


r(Xt)  = 


«>(*,•) 

w(X)’ 


If  the  integral  on  the  first  subinterval  is  narrower  than  its  share  of  the  total  tolerance 
requires,  then  the  tolerances  on  the  other  subintervals  are  relaxed  so  that  the  total  tolerance 


>Y-V‘V*' •• 


••  ■,  >  :  >  •• 
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can  be  met  more  efficiently.  On  the  other  hand,  if  the  integral  on  the  first  subinterval 
is  a  little  too  wide,  then  the  integrals  on  the  remaining  subintervals  can  sometimes  be 
commuted  accurately  enough  that  the  total  tolerance  can  still  be  satisfied. 

7.  Implementation  Details 

It  follows  from  the  above  discussion  that  realization  of  adaptive,  self-validating  quadrature 
routines  on  a  computer  requires  that  the  following  features  be  supported: 

(i)  Interval  computation  (arithmetic  operations  and  standard  library  functions). 

(ii)  Subroutines  for  automatic  generation  of  Taylor  coefficients. 

(iii)  Accurate  scalar  product  of  interval  vectors  (to  minimize  width  due  to  roundoff 
error). 

The  program  INTE  [9]  was  written  in  FORTRAN  for  the  Sperry  1100,  on  which 
only  (i)  [25]  and  (ii)  were  supported.  As  mentioned  above,  INTE  performs  self-validating 
numerical  integration  without  adaptation.  The  microcomputer  language  Pascal-SC  [22] 
and  the  ACRITH  package  for  the  IBM  370  series  of  computers  support  (i)  and  (iii),  to  which 
the  routines  (ii)  have  been  added.  The  implementation  described  in  |6]  for  Pascal-SC  is 
immediately  adaptable  to  ACRITH.  The  program  described  here  is  written  in  FORTRAN 
for  use  with  ACRITH.  In  particular,  the  following  operators  and  library  functions  are 
supported: 


+ 

SIN 

S1NH 

EXP 

- 

COS 

COSH 

LOG=LN 

* 

TAN 

TANH 

LOGlO 

/ 

COTAN=COT 

COTANH=COTH 

ERF 

**  constant 

ASIN 

ASINH 

ERFC 

ABS 

ACOS 

ACOSH 

GAMMA 

SQR 

ATAN 

ATANH 

LGAMMA 

SQRT 

ACOTAN=ACOT 

ACOTNH=ACOTH 

The  Taylor  series  terms  Fn(C,H)  and  Fn(X,H)  are  calculated  by  using  the  well- 
known  recurrence  relations  [15J,  |16j,  [19j  for  the  operators  and  functions  given  above.  In 
this  application,  the  “point”  of  expansion  and  the  stepsize  are  interval-valued  to  give  the 
desired  inclusions,  but  the  recurrence  relations  remain  the  same.  It  is  possible  for  the  user 
to  augment  the  list  of  library  functions  given  above  if  the  required  recurrence  relation  for 
Taylor  coefficients  of  the  new  function  is  known. 

Given  an  integrand  of  the  type  considered,  for  example, 

a.)  /<*>  =  rh' 

the  program  first  parses  it  into  a  code  list  [19]: 


Operator 

Operand  1 

Operand  2 

Result 

SQR 

x  (=Templ) 

Temp4 

4- 

1  (=Temp2) 

Temp4 

Temp5 

/ 

4  (=Temp3) 

Temp5 

F 

In  order  to  compute  the  series  for  /  expanded  at  C  with  stepsize  H,  this  code  list  is 
interpreted  to  obtain  the  sequence  of  calls 

Tempi  :=  (C,H)  {The  series  for  x.} 

Temp2  :=  (1)  {Constant  series.) 

Temp3  :=  (4)  {Constant  series.) 

Call  ITSQR(Templ,Temp4) 

Call  ITADD(Temp2.Temp4,Temp5) 

Call  ITDIV(Temp3,Temp5,F). 

The  result  is  an  array  which  contains  the  interval-valued  series  for  /  and  an  indication  of 
how  many  terms  were  computed.  For  example,  the  subroutine  ITSQR(t/,  V)  computes  the 
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series  for  V  =  U2,  given  the  series  for  U ,  by  means  of  the  recurrence  relations 


(7.2) 


i  -  l 

2  *  YjJj  *  Ui-j+x  +  lSQR(l/iii), 

*'/  2 

2  *  ^  ''Uj  *  Ui-j+i, 


i  odd, 


t  even, 


V  i  =  i 

where  V,  =  V,(C,  //)  and  l/t  =  Ui(C,H)  denote  the  tth  Taylor  coefficents  of  U  and  V, 
respectively  [19],  p.  49.  The  interval  function  ISQR  computes  ISQR(A')  =  X 2  instead  of 
X  ■  X ,  which  is  preferable  in  general,  since,  for  example,  {— 1,  l]2  =  [0, 1)  while  [-1,1]  • 
[— 1, 1]  =  [-1,1].  The  parsing  and  interpretation  at  runtime  described  above  is  unnecessary 
in  Pascal-SC,  because  the  compiler  generates  the  required  code  [6]. 

The  interval  Taylor  coefficients  F\(X,  C ),  F2(X,  C),  ...  are  maintained  in  a  record-like 
structure.  If  “F”  is  the  name  of  the  function  being  expanded,  then 
LF  Index  of  last  known  nonzero  term; 

MF  Index  of  last  known  term; 

OFL  Vector  of  series  terms — left  (lower)  bound; 

OFR  Vector  of  series  terms — right  (upper)  bound. 


The  designations  for  other  variables  replace  “F”  in  the  above. 

The  algorithms  used  depend  on  whether  the  endpoints  of  the  interval  X  of  integration 
are  elements  of  S,  that  is,  machine  numbers,  or  whether  they  are  intervals.  The  basic 
algorithm  is  for  the  case  X  =  ja,6]  €  IS,  and  is  called  the  RR  (REAL-REAL)  algorithm: 

1.  Compute  the  integral  on  [a,  6]; 

2.  Add  [a,  6]  to  the  list  of  subintervals; 

3.  Loop 

4.  Find  the  subinterval  on  which  the  width  of  the  integral  is  largest; 

5.  Bisect  it; 

6.  Compute  the  integral  on  the  left  subinterval; 


7.  Add  the  left  subinterval  to  the  list; 

8.  Compute  the  integral  on  the  right  subintervai; 

9.  Add  the  right  subinterval  to  the  list: 

10.  Compute  the  integral  on  [a,  6]  by  summing  the  integrals  on  all  the  subintervals; 

11.  Exit  when  accuracy  tolerance  is  met; 

12.  Exit  with  warning  when 

13.  no  further  improvement  in  accuracy  is  possible, 

14.  or  M  function  evaluations  are  exceeded; 

15.  End  loop. 

Each  subinterval  X  is  maintained  in  a  data  structure  of  the  following  form: 

XA  Left  endpoint  of  the  subinterval; 

XB  Right  endpoint  of  the  subinterval; 

OPTORD  Order  of  the  derivative  used  to  compute  the  remainder; 

WIDINT  Width  of  the  integral  on  this  subinterval; 

SINT  Interval-valued  integral  on  this  subinterval; 

WEGHT  Vector  of  interval  valued 

weights  (Gauss  and  Newton-Cotes), 
stepsize  (Taylor); 

FNVAL  Vector  of  interval-valued 

function  values  (Gauss  and  Newton-Cotes), 
series  terms  (Taylor); 

FNTRN  Vector  of  interval-valued 

function  values  (Gauss  and  Newton-Cotes), 
series  terms  (Taylor), 
including  remainder  terms. 

At  steps  1,  6,  and  8,  the  integral  is  computed  on  the  subinterval  [XA,XB]  using  one- 
panel  versions  of  Gauss,  Newton-Cotes,  or  Taylor  polynomials  as  outlined  in  Sections  2  and 
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3.  The  weight  vectors  and  functions  are  arranged  in  the  vectors  WEGHT  and  FNTRN, 
respectively,  in  such  a  way  that  the  interval  inclusion 

(7.3)  J  =  WEGHT,  *  FNTRN, 

i 

of  Sa  f(X)dl  is  computed  as  a  single  inner  product.  Similarly, 

(7.4)  Rf  =  WEGHT*  *  FN VALt. 

t 

At  step  13,  if  J  C  Rf,  then  the  loop  is  exited.  In  this  case,  further  reduction  of  the  width 
of  the  truncation  error  cannot  reduce  w(J).  For  Gauss  and  Newton-Cotes  formulas,  SINT 
is  used  only  to  give  WIDINT.  For  Taylor  polynomials,  SINT,  is  intersected  with  (7.3). 
SINT,  is  computed  using  the  intersection  principle  discussed  at  the  end  of  §3.  For  a  few 
subintervals,  ^SINT,  is  narrower  than  (7.3),  while  the  situation  is  usually  reversed  for  a 
large  number  of  subintervals,  because  (7.3)  uses  the  accurate  scalar  product. 

The  arrangement  of  subintervals  in  the  arrays  listed  above  must  be  relatively  straight¬ 
forward  to  allow  J  to  be  computed  by  a  single  scalar  product  operation.  Each  iteration  of 
the  loop  from  step  3  to  step  15  removes  one  subinterval  from  the  list  and  replaces  it  with 
two  subintervals.  Subintervals  are  not  otherwise  deleted  from  the  list.  Hence,  the  follow¬ 
ing  simple  allocation  scheme  works:  On  the  tth  pass  through  the  loop,  the  information 
about  the  left  subinterval  is  stored  in  the  locations  previously  used  by  its  parent,  and  the 
information  about  the  right  subinterval  is  stored  in  the  (t  +  l)st  locations,  following  the 
already  computed  values.  Hence,  insertion  requires  no  searching.  The  widest  subinterval 
is  found  at  step  4  by  a  sequential  search  of  the  array  WIDINT(l..t). 

By  contrast,  QUADPACK  [17]  maintains  its  list  of  pending  subintervals  in  sorted 
order,  so  no  search  is  necessary  for  the  next  subinterval  to  be  processed.  However,  new 
subintervals  are  inserted  at  locations  found  by  a  sequential  search,  followed  by  changing 
pointers  to  all  following  entries  in  the  list.  For  each  subinterval  processed,  the  program 
does  one  sequential  search,  while  QUADPACK  does  two.  In  addition,  QUADPACK  uses 
two  sets  of  pointer  adjustments. 
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For  integration  using  Taylor  polynomials,  the  maintenance  of  list  of  subintervals  is 
somewhat  more  complicated,  because  the  program  reuses  the  series  which  it  has  previ¬ 
ously  computed.  To  illustrate  the  ideas,  consider  the  first  execution  of  the  loop  at  step  3. 
At  that  point,  the  list  of  subintervals  contains  only  one:  X  =  [a,  6)  itself.  FNTRN  contains 
OPTORD-1  terms  of  the  series  for  /  expanded  at  c  =  (a  +  b)j 2  and  the  truncation  error 
term  involving  /tIoptord)^)  gjven  by  equation  (4.6).  Provided  that  the  requested  tol¬ 
erance  exceeds  the  noise  inherent  in  the  function  evaluation,  a  stepsize  h  can  be  computed 
which  is  small  enough  that  the  requested  tolerance  per  unit  step  is  satisfied  on  the  interval 
[c  -  h,c  -r  A].  Notice  that  if  a  relative  tolerance  is  requested,  then  this  requires  a  current 
estimate  for  J .  The  value  of  the  integral  on  this  subinterval  can  be  computed  at  a  cost 
proportional  to  OPTORD,  instead  of  a  cost  proportional  to  OPTORD2,  which  would  be 
required  to  generate  J  directly  by  using  the  recurrence  relations  for  /.  Following  this,  the 
two  subintervals  |a,c  -  A]  and  [c  +  A,  6]  are  processed  directly.  Consequently,  this  method 
breaks  the  subinterval  of  integration  into  three  parts,  rather  than  bisecting  it. 

Thus,  for  integration  by  Taylor  polynomials,  step  5  of  the  RR  algorithm  is  replaced 

t>y 

5.0'  Compute  h  such  that  the  tolerance  is  satisfied  on  [c  -  h,c  -f  h J; 

5.1'  Compute  the  integral  on  \c  -  h,c  +  h]  from  information  in  FNVAL; 

5.2'  “left  subintervar  :=  [XA.  c  -  h\\ 

5.3'  “right  subinterval”  :=  jc  +  h,  XBj. 

The  middle  subinterval  jc  -  h,c  +  h]  is  maintained  on  the  list  of  subintervals  so  that 
ts  contribution  to  J  is  included  in  the  scalar  product  in  (7.3).  This  has  the  helpful  side 
ffect  that  h  can  be  chosen  somewhat  optimistically.  If  the  choice  is  too  optimistic,  then 
:  -  h,c  +  h\  will  be  selected  later  for  further  processing  as  the  worst  subinterval,  at  which 
ime  it  will  be  broken  into  three  parts.  One  of  the  new  subintervals  will  occupy  the  place 
f  the  parent  interval  in  the  list  maintained,  while  the  other  two  will  be  added  to  the  end 
f  the  list. 
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A  further  refinement  could  be  implemented.  The  stepsize  h  computed  at  step  5.0'  has 
the  following  property:  On  each  subinterval  of  length  2 h  which  is  contained  in  [XA,XBj, 
the  use  of  a  Taylor  polynomial  of  degree  OPTORD-1  yields  an  integral  which  satisfies  the 
requested  tolerance.  The  integration  on  such  a  subinterval  can  be  done  with  half  the  usual 
work.  No  series  for  the  truncation  error  needs  to  be  computed  because  the  truncation 
error  can  be  bounded  by  using  the  global  remainder  term  on  [XA,XB].  If  [XA,XB]  can 
be  covered  by  a  few  subintervals  of  length  2 h,  then  this  could  be  done,  and  division 
into  three  parts  would  be  needed  only  when  the  middle  part  is  relatively  small.  This 
refinement  could  improve  the  efficiency  of  the  program.  However,  the  improvement  would 
likely  be  modest,  because  the  advantages  of  intersecting  subsequent  estimates  on  each 
small  subinterval  would  be  lost,  and  the  number  of  subintervals  added  to  the  list  would 
no  longer  be  constant.  The  program  also  does  not  reuse  function  evaluations  required  by 
Gauss  or  Newton-Cotes  formulas,  although  it  could  be  modified  to  do  so. 


8.  Numerical  Examples 

We  give  four  examples  to  illustrate  the  accuracy  and  the  reliability  of  the  program.  All 
computations  were  done  in  double  precision  on  an  IBM  4341  computer,  using  calls  to 
ACRITH  routines  for  all  necessary  interval  calculations,  including  scalar  products. 

[0.7  j 

Example  1.7=  /  - dx. 

Jo. C  1  -  x 

Interval  bounds  for  the  answer  can  be  computed  in  three  ways: 

(8.1)  7  =  log(l  —  0.6)  -  log(l  -  0.7), 


(8.2)  7  =  log(4/3), 

or  by  adaptive,  self-validating  quadrature. 

Neither  0.6  nor  0.7  are  machine  numbers,  so  they  are  converted  to  intervals  which 
are  one  machine  number  wide.  All  three  methods  give  the  interval  |0. 2876820724517808. 
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876820724517811],  but  the  results  are  20,  13,  and  14  machine  numbers  wide,  respec- 
;ly.  That  is,  the  program  is  capable  of  accuracies  comparable  to  evaluation  of  the 
dytic  expression  for  the  answer.  This  result  required  43  equivalent  function  evaluations 
18  subintervals.  In  order  to  appreciate  the  accuracy  achieved  by  the  program,  we  must 
sider  some  details  at  the  level  of  machine  numbers.  If  (8.1)  is  evaluated  without  sim- 
ication  using  interval  calculations,  the  result  is  the  hexadecimal  interval  [Z  4049  A588 
>3  6E41,  Z  4049  A588  44D3  6E55j,  which  is  20  machine  numbers  wide.  Using  (8.2),  the 
4  hexadecimal  digits  are  jZ  ...  6E45,  Z  ...  6E52],  which  is  13  machine  numbers  wide. 
>  adaptive,  self-validating  quadrature  program  gives  jZ  ...  6E44,  Z  ...  6E52],  which  is 
nachine  numbers  wide. 


Example  2. 


(see  §5.1). 


This  example  illustrates  the  order  adaptation  required  to  handle  the  nonexistance  of 
).  The  program  gave  the  interval  [5.33333  33333  27,  5.33333  33333  36],  It  stopped  when 
id  used  400  effective  function  evaluations  on  226  subintervals.  Most  of  the  subintervais 
!  clustered  near  the  origin.  Away  from  0,  as  many  as  13  series  terms  were  used. 


Example  3. 


s: 


f(x)dx,  where  f(x) 


0,  i  <  0.3, 
1  i  >  0.3. 


This  example  illustrates  integration  of  a  simple  discontinuous  function.  The  parser 
not  accept  a  piece-wise  definition,  so  this  function  was  coded  by  hand.  Table  1  shows 
esults  for  tolerances  0.0  and  1.0E-15,  and  for  recognizing  that  the  function  has  a 
larity  on  the  interval  of  integration. 
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Width  of 

Integral 

Function 

Evaluations 

Subintervals 

Tolerance  =  0.0 

|0, 1] 

2.50E-16 

218 

146 

(0,  0.3]  +  |0.3,  1] 

2.78E-17 

22 

29 

Tolerance  =  1.0E-15 

|0,1| 

7.77E-16 

194 

130 

|0,  0.3]  +  [0.3,  1] 

1.39E-16 

10 

5 

Table  1.  Integrating  a  Discontinuous  Function. 


As  is  usually  the  case,  the  program  performs  much  better  when  the  user  recognizes 
the  presence  of  a  discontinuity.  The  performance  would  be  even  better  if  the  discontinuity 
occured  at  a  machine  number.  Notice  that  the  cost  of  requesting  the  program  to  do  the 
best  it  can  (tolerance  =  0.0)  is  only  slightly  more  than  the  cost  when  a  lesser  tolerance 
is  prescribed.  The  table  shows  more  subintervals  than  evaluations  because  evaluating  a 
series  which  is  =  0  is  not  counted  as  an  evaluation. 

fl  1 

Example  4.  /  - r dx. 

Jo  1  -  ox2 

This  example  illustrates  the  performance  of  the  program  as  the  integrand  varies  from 
very  smooth  to  being  undefined  at  a  point  in  the  interval  of  integration.  These  calculations 
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done  in  single  precision  with  an  absolute  error  tolerance  request  of  1.0E-5. 


Error 

Function 

Maximum 

Absolute 

Alph 

a  Code 

Evaluations 

Subintervals 

Order 

Error 

0.00 

0 

2 

2 

3 

0.0 

0.05 

0 

8 

2 

20 

4.4E-16 

0.10 

0 

8 

2 

20 

1.8E-14 

0.15 

0 

8 

2 

20 

6.6E-12 

0.20 

0 

8 

2 

20 

7.8E-10 

0.25 

0 

5 

2 

15 

4.1E-06 

0.30 

0 

15 

6 

15 

1.9E-10 

0.35 

0 

12 

6 

13 

4.5E-08 

0.40 

0 

12 

6 

13 

1.8E-07 

0.45 

0 

21 

10 

15 

6.1E-08 

0.50 

0 

19 

10 

14 

1.4E-07 

0.55 

0 

16 

10 

12 

7.9E-06 

0.60 

0 

24 

14 

12 

2.4E-07 

0.65 

0 

35 

18 

15 

6.4E-08 

0.70 

0 

29 

18 

12 

5.0E-06 

0.75 

0 

38 

22 

12 

1.3E-07 

0.80 

0 

48 

26 

15 

6.6E-08 

0.85 

0 

43 

26 

12 

7.5E-06 

0.90 

0 

62 

34 

15 

6.8E-08 

0.95 

0 

73 

42 

14 

1.2E-07 

1.00- 

66 

254 

150 

8 

3.2E-01 

1.05 

67 

2 

Table  2. 

Performance  Profile. 

error 

code  of  66 

signals  that  the  program  was 

unable  to  meet 

the  requested  tol- 

while 

67  means 

that  it  was  unable  to  evaluate  the  integrand 

As  the  problem 

ne  more  difficult,  the  program  required  more  effective  function  evaluations  and  more 
tervals. 
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